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METHODS AND SYSTEMS FOR MODEL REDUCTION AND SYSTEM 
IDENTIFICATION OF DYNAMIC SYSTEMS WITH MULTIPLE INPUTS 

Inventor 
Taehyoun Kim 

10 



Field of the Invention 

This invention relates generally to methods and systems for efficient model reduction 
15 and system identification of dynamic systems, and more specifically, for linear dynamic 
systems having multiple inputs using a new Eigensystem Realization Algorithm. 



Background of the Invention 

20 Most modern dynamic models are constructed based on finite spatial discretizations 

of continuous systems, often resulting in a considerable number of degrees of freedom in the 
model. Consequently, for fast and efficient estimation of dynamic behavior, as well as 
optimizations and closed-loop control designs, , a model reduction is typically performed. 
Desirable characteristics of a reduced-order dynamic model include that the size of the 

25 system must not be too large, the model must be robust and have a good fidelity, it must be in 
the state-space, time domain formulation for implementation of active control systems and 
nonlinear time analysis, and finally, the reduction process itself must not be too expensive. 

There have been many model reduction methods available, but most of them require 
modifying the main frame of the computational model, and are prone to a long model 
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construction time if the model is subjected to many driving inputs. The latter in particular is 
true in unsteady Computational Fluid Dynamic (CFD) applications where the moving solid 
boundary is often described by many structural mode inputs. For example, a typical 
commercial airplane simulation may involve as many as 200 structural modes. 
5 Recently, Eigensystem Realization Algorithm (ERA) was used successfully in 

aeroelastic flutter predictions using a computational fluid dynamics code (CFD) developed 
by the NASA Langley Research Center known as CFL3D. (See Juang, J.N., Applied System 
Identification, Prentice Hall Englewood Cliffs, New Jersey, 1994; Silva, W.A., and Bartels, 
R.E., Development of Reduced-Order Models for Aeroelastic Analysis and Flutter Prediction 

10 Using CFL3Dv6.0 Code, AIAA-2002-1596; and Hong, M.S., Bhatia, K.G., SenGupta, G., 
Kim, T., Kuruvila, G., Silva, W.A., Bartels, and R., Biedron, R., Simulations of a Twin- 
Engine Transport Flutter Model In the Transonic Dynamics Tunnel, IFASD Paper 2003-US- 
44). The ERA method, which is usually used as a system identification technique, has a very 
attractive feature in that unlike model reduction methods based on Galerkin scheme (e.g., 

15 Dowell, E.H., Hall, K.C., and Romanowski, M.C., Eigenmode Analysis in Unsteady 
Aerodynamics: Reduced-Order Models, Applied Mechanics Review, Vol. 50, No. 6, 1997, 
pp. 371-386), there is no need for on-line implementation of the algorithm. That is, the ERA 
method is a post-processing tool that identifies and generates system matrices based on the 
input and output data alone. 

20 FIGURE 1 is a schematic view of the ERA (a.k.a. Pulse/ERA) method 100 of 

modeling a dynamic system in accordance with the prior art. For simplicity, only its 
fundamental state-space realization theory which is attributed to Ho and Kalman is discussed, 
(see Ho, B.L. and Kalman, R.E., Effective Construction of Linear State-Variable Models 
from Input/Output Data, Proceedings of the 3rd Annual Allerton Conference on Circuit and 

25 System Theory, 1965, pp. 449-459). For a general description of the Pulse/ERA method, see 
Juang, J.-N. and Pappa, R.S., An. Eigensystem Realization Algorithm for Modal Parameter 
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Identification and Model Reduction, Journal of Guidance, Control, and Dynamics, Vol. 8, 
1985, pp. 620-627. 

In the Pulse/ERA method 100, it is assumed that the dynamic system under 
consideration can be described in the following finite-dimensional, discrete-time form: 



wibere 



x s (X x 1} state TOctar 
lai = (jYi x 1) input rector 
y s (jV ff x i) output weter 



(1) 
(2) 



(3) 
(4) 
(5) 



The system matrices (A, B), (A, C) are assumed controllable and observable. First, given M 

+ 1 equally distributed time steps, s 3:1 ^ ( tt = 0 b 1 i 2 i-"bAO j( for a single i-rt input 
10 vector the system output is sampled subjected to the unit pulse, 



j = p m I 1 (a = 0) 1 



(6) 



Assuming zero initial condition, x° = x(0) = 0, one obtains 
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yf = CA jlM l^ (T) 



5, I-th cxflmiuiof B (S) 
s i-th cdhuim of D (9) 

The constant matrices in the above sequence are known as the Markov parameters. As shown 
in FIGURE 1, for all inputs 102 into the dynamic system, the Markov parameters are 
5 computed (block 104) creating the sequence of pulse-response matrices: 

v = W ... y%\ (» -otW-JO (io) 

Next, based on the system Markov parameters (block 104), the Pulse/ERA method 100 
defines N 0 x (Ni x(M — 1)) Hankel matrices: 

Hfl = fY l Y 2 ... Y^'H 

bb ClJ A A 2 ... A^B (11) 
Hi = fY 2 Y 3 ... Y^\\ 

- CjjA At ... A M - l |jB (12) 

10 

Singular-value-decomposition (SVD) of H 0 yields 
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- U^-sJfvS (13) 



where R = rank(H 0 ). Finally, a balanced realization of the system under question is obtained 
by pseudo-inverting various submatrix components (block 108 of FIGURE 1) as follows: 



D = Y D (%xl) (14) 

C ^ U*s)f (jV.xJQ (15) 

B ~ the first jY; odiums of Sjj- {R x jY*) (16) 

A * S^-UlH L V fi S^ fJ£xJfi (ID 



Since R « L, the above model represents a reduced-order realization of the original 
system. Note that the realization is optimal in that it is balanced between inputs and outputs. 
10 However, the total number of samples taken is Nt x (M+ 1), which increases proportional to 
the number of inputs. Also, for an accurate system identification Ho must have sufficient 
columns and rows, i.e., Nt x (M- 1), >R and N 0 >/?. 

For a very large data set Y n with many time steps and a large number of inputs, the 
Eigensystem Realization Algorithm/Data Correlations (ERA/DC) method may be preferred. 
15 As described, for example, in Juang, J.N., Applied System Identification, Prentice Hall 
Englewood Cliffs, New Jersey, 1994, the ERA/DC method may be used to compress the 
amount of data and reduce the computation time required for the SVD of the Hankel matrix. 

Although desirable results have been achieved using the prior art computational 
methods, there is room for improvement. For example, if the unsteady CFD model is driven 
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by multiple structural inputs, as described above, the computation time required to obtain all 
the pulse responses increases proportional to the number of the inputs, making the ERA 
method undesirably slow and inefficient. Therefore, a need exists for improved methods for 
model reduction and system identification of large-scaled linear dynamic systems having 
5 multiple inputs 



Summary of the Invention 

The present invention is directed to methods and systems for model reduction and 

10 system identification of dynamic systems having multiple inputs. Embodiments of the 
methods in accordance with the present invention may advantageously reduce the model 
construction and system identification time, especially for a large-scaled system, and 
compress the amount of output data. Furthermore, they do not require modifying the original 
code, and take only input and output data for the model construction. That is, it is 100% a 

1 5 post-processing tool. 

In one embodiment, a method of model reduction and system identification includes 
generating a plurality of statistically independent random numbers for use as input signals, 
and performing a singular-value-decomposition directly on the system response due to the 
simultaneous excitation of all the inputs, also known as the Single-Composite-Input (SCI). 

20 In alternate embodiments, the method further includes sampling individual pulse responses 
for the first time and second time steps. Alternately, the input signals are filtered through a 
low-pass filter. The plurality of input signals may further include applying multiple step 
inputs in a sequential manner, and applying multiple pulse inputs in a sequential manner. 

In an alternate embodiment, a method of model reduction and system identification 

25 includes generating a plurality of statistically independent random numbers for use as input 
signals, performing a singular-value-decomposition directly on the system response due to 
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the Single-Composite-Input (SCI), sampling individual pulse responses for the first time and 
second time steps. The method further includes defining the Hankel-like matrices Hco and 
H c i as follows: 



10 



= Cfc x 3tT ... (25) 

Hd s bii & ... 

- CAf«* 3T ... (26) 

= UjtSSf S# 3 vE (27) 

Finally, the method identifies the system matrices (A, B, C, D) by a least-square 
approximation as follows: 

D - Y* (28) 

C a U^S^ 2 (29) 

B a E^U^Y 1 (30) 

A c= E^I^H^VjiE^ (31) 



Brief Description of the Drawings 

The preferred and alternative embodiments of the present invention are described in 
detail below with reference to the following drawings. 
15 FIGURE 1 is a schematic view of the Pulse/ERA method in accordance with the prior 

art; 
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FIGURE 2 is a schematic view of the SCI/ERA method in accordance with an 
embodiment of the present invention; 

FIGURE 3 is a schematic view of a vortex lattice model for computing an unsteady 
flow field around a flat, rectangular wing in accordance with an embodiment of the present 
5 invention; 

FIGURE 4 shows three sets of random numbers generated as inputs to the vortex 
lattice formulation of FIGURE 3; 

FIGURE 5 shows sets of generalized aerodynamic forces computed using the vortex 
lattice formulation of FIGURE 3; 
10 FIGURE 6 shows a set of aeroelastic eigenvalues computed using the vortex lattice 

formulation of FIGURE 3 at a first scale for V = 80 m/sec; 

FIGURE 7 shows a set of aeroelastic eigenvalues computed using the vortex lattice 
formulation of FIGURE 3 at a second scale for V = 80 m/sec; 

FIGURE 8 shows an aeroelastic time response of the first structural mode computed 
1 5 using the vortex lattice formulation of FIGURE 3 for V = 80 m/sec; 

FIGURE 9 shows an aeroelastic time response of the second structural mode 
computed using the vortex lattice formulation of FIGURE 3 for V = 80 m/sec; 

FIGURE 10 shows a set of aeroelastic eigenvalues computed using the vortex lattice 
formulation of FIGURE 3 at a first scale for V = 121.2 m/sec; 
20 FIGURE 1 1 shows a set of aeroelastic eigenvalues computed using the vortex lattice 

formulation of FIGURE 3 at a second scale for V = 121.2 m/sec; 

FIGURE 12 shows an aeroelastic time response of the first structural mode computed 
using the vortex lattice formulation of FIGURE 3 for V = 121.2 m/sec; 

FIGURE 13 shows an aeroelastic time response of the second structural mode 
25 computed using the vortex lattice formulation of FIGURE 3 for V = 121.2 m/sec; 
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FIGURE 14 shows a graph of model construction time of different ERA methods 
versus number of inputs for the vortex lattice formulation of FIGURE 3; 

FIGURE 15 shows a graph of CPU seconds spent constructing various ERA models 
in accordance with alternate embodiments of the present invention (on which FIGURE 14 is 
5 based); 

FIGURE 16 is a representative CFL3D simulation model in accordance with an 
embodiment of the present invention; 

FIGURE 17 shows four structural mode shapes used for the CFL3D simulation model 
of FIGURE 16; and 

10 FIGURE 18 shows a V-g plot of two different reduced-order aeroelastic models based 

on the CFL3D simulation in accordance with an embodiment of the present invention. 



Detailed Description 

15 The present invention relates to methods and systems for model reduction and system 

identification of dynamic systems having multiple inputs. Many specific details of certain 
embodiments of the invention are set forth in the following description and in FIGURES 2- 
18 to provide a thorough understanding of such embodiments. One skilled in the art, 
however, will understand that the present invention may have additional embodiments, or 

20 that the present invention may be practiced without several of the details described in the 
following description. 

In general, embodiments of methods in accordance with the present invention include 
a new, efficient discrete-time domain system identification and model reduction method for 
large-scaled linear dynamic systems with multiple inputs. In one embodiment, a method in 

25 accordance with the present invention is based on a modification of the classical Eigensystem 
Realization Algorithm (ERA) and a simultaneous injection of multiple inputs, so called the 
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Single-Composite-Input (SCI). Since the system response is sampled almost exclusively for 
the single representative input instead of the multiple individual inputs, embodiments of the 
present invention can significantly reduce the model construction time as well as the amount 
of the sampled data. Derivation of the new algorithm in accordance with the present 
5 disclosure may include performing a singular value decomposition using a single set of 
output measurements that are not necessarily attributed to pulse inputs. Representative 
simulations performed using embodiments of the present invention exhibit reduced 
computation times and excellent results obtained from the reduced-order models, thereby 
showing great potential of the present invention as a linear system identification and model 
10 reduction tool for large-scaled systems subjected to multiple inputs. 

Unless otherwise stated, the following nomenclature is used throughout the following 
detailed description: ' 



A B C D System matrices 

1 5 Adi A d 2 Aeroelastic system matrices 

A s B s C s Structural system matrices 

b Reference length 

Qj Cross-correlation coefficient 

E Expected value 

20 K Covariance matrix defined in (48) 

. L Dimension of original system 

M Number of time or frequency samples 
or Mach number 

m c k Mass, damping, stiffness matrices 

25 p (Rjx 1) generalized coordinate vector 

q Dynamic pressure (=1/2/? V 2 ) 

R Number of chosen singular modes or 
the dimension of realized model 
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Ri Number of chosen KL modes 

s Laplace variable 

t Real time 

u (Nixl) input or generalized structural 

5 coordinate vector 

Ui i-th input or structural coordinate 

V Air speed 

x (L x 1) state or aerodynamic state vector 

X Fourier transform of x 

1 0 y (N 0 xl) output vector or (NjXl) 

generalzied aerodynamic force vector 

y n i Pulse response due to i - th input 

Oi KLmode 

* KL modal matrix 

15 p Air density 

x Reduced time ^ 

© Frequency (rad/sec) 

co c Maximum cut-off frequency 

20 Subscripts 

i Input 

o Output 

R Reduced 

ref Reference 

25 s Structure 



FIGURE 2 is a schematic view of the Single-Composite-Input (SCI)/ERA method 
200 in accordance with an embodiment of the present invention. As described more fully 
below, the SCI/ERA method 200 is based on the premise that for a linear system, one can 
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apply the multiple inputs simultaneously and get all the system responses that are necessary 
for the model reduction. Since the computational model needs to be executed almost 
exclusively for the representative input, the model construction time can be significantly 
reduced. 

Embodiments of the present invention include implementing the SCI/ERA method for 
fast and efficient model reduction of linear, finite-dimensional, discrete-time systems with 
multiple inputs. To accommodate the SCI within the framework of the Pulse/ERA, it is 
necessary to modify the original algorithm. In particular, as shown below, the new 
formulation does not rely on the system Markov parameters explicitly. Instead, it performs 
the singular-value-decomposition (SVD) directly on the output measurements that are in 
general not attributed to pulse inputs. Statistically independent random numbers are used in 
lieu of the pulses for the multiple input signals. Naturally, embodiments of the present 
invention can also be used towards system identification provided that all the time 
measurements are available from experiments. Application of the SCI/ERA method 200 to 
computational fluid dynamic systems and formulation of reduced-order aeroelastic models 
are presented below, where it is shown that depending on how the displacement and velocity 
inputs on the moving boundary are treated, two different kinds of reduced-order aerodynamic 
and aeroelastic models may be generated. 

For a demonstration of the proposed method, two CFD models are considered in the 
examples below. They are the Vortex Lattice Model (VLM) for inviscid, subsonic, 
incompressible flow, and the CFL3D for viscous, transonic, compressible flow. Reduced- 
order aeroelastic models are also constructed by combining the reduced aerodynamic models 
and the structural system. It is shown that not only does the new method shorten the model 
construction time substantially, but the accuracy of the resulting reduced-order models is 
excellent. The proposed scheme has a great potential as a linear system identification and 
model reduction tool for large-scaled systems subjected to multiple inputs. 
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With reference to FIGURE 2, the SCI/ERA method 200 proceeds as follows. First, at 
a block 202, individual pulse responses are sampled for the first two time steps: 



Y D - M y§ ... 4-j (18) 
Y l = hi y\ ... yy (19) 

Next, at a block 204, an SCI is constructed as follows: 



wluens 



tf?<7 s I>*f (for states) (20) 
' d&y = £darf (for outputs] (21) 



j£ = a secLueuce of arbataoy numbers (22) 



10 To ensure independency of the inputs, it is desirable that the signals be as 

uncorrected as possible. In an ideal case they would be statistically uncorrected random 
signals, i.e., Q/m) =E[r n x r Y7 = 0 for i but they are hard to construct for numerical 
analysis. 

Subject to the SCI, at a block 206, the SCI/ERA method 200 samples the system 
15 response y" for n = 0, 1, 2, . . . , M t to obtain: 



- y"-Erftf (23) 

= ^-rrf^-Ssifl (24) 
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Note that ^os^i are measurements of the states in the reduced dimension of C and C 

A. 

Next, at a block 208, similar to the Hankel matrices, the SCI/ERA method 200 
defines the following Hankel-like matrices: 

= cix 1 xr ... ^-'[i (25) 
Hd = hii yf 3 ... y&~H 

= CAlx 1 x 2 ... x itf - , i (26) 

STVD of H^ji yi^Ms 

He ■ USV 

= U*sf2$ 2 VE (27) 

where = rank(H C0 ). The size of the above matrices is NoX (M - 1), N, times smaller than 
the previous Ho, Hi defined in the Pulse/ERA. 

Thus, at a block 210, the new realization is then 



D = Y° (28) 
C * U*Sif (29) 
E cr SJ^U^Y 1 (30) 
A st Sj^USH^V*!^' 2 (31) 

Unlike the Pulse/ERA method 100 described above with reference to FIGURE 1, the 
SCI/ERA method 200 may be optimal in that it may be balanced between states and outputs. 
15 As in the previous method 100, for an accurate realization H co must have (M - 1) >R and N 0 
^R. However, using the SCI/ERA method 200 in accordance with the present invention, the 
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total number of samples taken may be only M + 1 +2 x which is much less than the 
previous Ni x (M+ 1) when M samples of pulse response are taken for each input. 

It will be appreciated that compared to the prior art Pulse/ERA method 100, 
embodiments of the present invention may advantageously require a much smaller set of time 
5 measurements, thereby reducing significantly both the computation time and the bulk of the 
output data. Furthermore, the new H C09 H ci are constructed based on the sampled system 
states subjected to combined random inputs, and as such they are not directly related to the 
Markov parameters. However, at block 202 embodiments of the present invention (e.g. the 
SCI/ERA method 200 shown in FIGURE 2) do require the first two pulse responses yf and 
10 yt for each input in order to estimate the state measurements according to Equations (23) and 
(24). 

Although the use of random signals is described, other types of signals can also be 
used for the SCI provided that they are statistically independent. For a linear system, any 
arbitrary response contains the fundamental characteristics under the assumption that the 

15 system is controllable and observable. Embodiments of the present invention may 
advantageously employ this principle, along with the principle of superposition. 

An alternate embodiment to the SCI/ERA method discussed above will now be 
described. As stated above, an important requirement in any ERA method is that to improve 
the accuracy of the model construction, one must have a sufficient number of output 

20 measurements, more than the number of singular modes that are extractable from the Hankel 
or Hankel-like matrices. Failure to satisfy this requirement implies that we don't have enough 
modes to approximate the system response. When the requirement is not met, assuming again 
that (A, C) is observable we can expand the data matrices by sampling additional responses 
as follows: 

25 
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Hen, = 



" c 






CA 






CA 5 

■ 


be 1 x 2 ... x-"- 1 ! 


■ 
■ 

CA* 






jcq jco — jcu 
JcX >cl ■"• Jel 








c 






CA 






CA 2 

■ 


A» JpC^* ■ ■ ■ ■ 




■ 

CA* 






til 


y&2 ■ • ■ yis 




J*ff4l • • • JcA'-lL J 



(32) 



(33) 



where 



& = CA h x" 
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A-hl jYj 



SVD of tie new H^u iksadis to 



= |Ua* U ai? j 



Sir 0 
0 0 



IK 1 



from vladii -we obtain 



D 
C 
B 
A 



= the first jV P rows of Yf 

c£ the first jY e rovra of Uj/jE}^ 



where 



yf 



(34) 



(35) 



(36) 
(37) 
(33) 
(39) 



(40) 



Using this alternate approach, the total number of measurements available is now (K 
+ 1) x No. It is noted that the additional time samples are required for the pulse response as 
well as for the response due to the SCI. More specifically, if K blocks of outputs are to be 
5 added pulse responses due to each input must be sampled at K additional time steps in 
addition the first two time steps, f = 0 audi t = At, Also, the response due to the SCI must 
be sampled at K additional steps beyond the M-th step. The total number of samples to be 
taken is thus M +1+ K + (2+ K) x N ( . K should desirably be sufficiently large enough to 
satisfy the measurement requirement, (K + 1) x N 0 Needless to say, this alternate 

10 scheme requires extra computation time because of the additional time samples required in 
the data set. 
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A few candidate signals for the SCI will now be described. For an ideal linear model 
any of these SCIs can be used and all of them should yield ROMs of essentially the same 
quality. 

As noted above, random signals may be used for construction of SCI. This approach 
5 may be referred to as random inputs SCI (RSCI). That is, use 



in Equations (20) and (21) above. 

Alternately, a filtered inputs based SCI (FSCI) may be employed. More specifically, 
10 to inject smooth inputs, one can filter the random signals through a low-pass filter. That is, 



Using such low frequency signals may allow better convergence when applying the 
SCI to CFD models. Furthermore, since the frequency content is limited, it is possible to 

15 generate a smaller ROM directly from the SCI/ERA without any further reduction. A 
potential drawback is that if the filtered signals become too narrowly banded, they may not 
be as uncorrelated as desired. However, based on the theory of ergodicity [9], the statistical 
independence could be fortified by using longer signals and sampling the response for a 
longer period of time, (see, e.g., Papoulis, A., Probability, Random Variables, and 

20 Stochastic Processes, McGraw-Hill Book Company, New York, New York, 1982, 
incorporated herein by reference). 

In another approach, a step input based SCI (SCCI) may be used. In analogy to a 
single step input, one can apply multiple step inputs in a sequential manner: 



a aequeiupe of random numbers 



(41) 




s a sequeiios of filtered! -random anunbeis 



(42) 
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If = 



= a step Impart ajjpEed at k^ih step (43) 

To assure independence of the inputs, the onsets of the signals should be apart from 
each other by a sufficient number of steps, i.e., k{ must be large enough. 
5 Similarly, in another alternate approach, a pulse input based SCI (PSCI) may be used. 

In this approach, one can also apply multiple pulse inputs in a sequential manner: 



= a step 'input appSed at fe-tli step (45) 

where 



10 Again, ^ should be large enough to ensure independency of the applied signals. 

Yet another alternate embodiment to the SCI/ERA methods discussed above will now 
be described. In applications of discrete-time computational models, there exist two 
conflicting requirements for the incremental time step At. For numerical convergence, one 
should adopt a sufficiently small At\. On the other hand, given the highest frequency of 

15 interest, coc, the Nyquist criterion requires Usually, for practical purposes 

Atf- » &Z\. 
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For instance, in a structural model that is coupled with a CFD model for aeroelastic 
applications, the highest mode usually has a much lower natural frequency than the highest 
frequency content in the aerodynamic model. If the signals used in the SCI/ERA methods are 
sharp as in the random, step, or pulse inputs, the SCI will excite all the system dynamics and 
5 hence this characteristic will be carried over to the ERA based reduced-order model As a 
result, the ERA reduced-order model (ROM) may be prone to have a large dimension to span 
the same high frequency range as the original full-order model (FOM), which suggests that 
there may be room for further order reduction in the ROM. 

To perform a second order reduction, one can apply the Frequency-Domain 

10 Karhunen-Loeve (FDKL) method to the ERA ROM defined by matrices, (28)-(31), wherein 
frequency samples of the system within the given frequency range, (-co 0 co c ) are used to 
extract optimal modes, and a new reduced-order model is constructed via the Galerkin's 
approximation, (see, e.g., Kim, T., Discrete-Time Eigen Analysis and Optimal Model 
Reduction for Flutter & Aeroelastic Damping/Frequency Prediction Based On CFL3D-ERA, 

15 B- AD VTECH-LLL-M02-0 1 3 , BCAG, February 2003, incorporated herein by reference). In 
this embodiment, the optimal or so called KL modes 0j, are the eigenmodes of the 
covariance matrix K: 



- As* (47) 

where 

Kij a XU)X(u^ r (48) 
= sainpfing frequences 

= |wl as uuujj (49) 
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where <£>\ = -co c and o M = co c . X(cc>i) are the frequency solutions of the ERA ROM 
subjected to the single-composite-input described by (20) and (21) except that it is prescribed 
in the frequency domain. Once the optimal modes are obtained, assume 



x ^ *p (150) 

where 

Pi 



J* 
■ 



(52) 



/?/ is set to be equal to the rank of the covariance matrix which is usually smaller than 



R. 



After inserting (50) into (1) and (2) with the ERA ROM matrices, premultiplying both 

T 

1 0 sides by <I> yields a new reduced-order model as 

p"" 1 ' 1 = Aip* + Bad" (53) 
y" = Cap" + 0111" (54) 

Aa s * T A* (55) 
Bj = * T B (56) 
C, = C-$ (57) 

The dimension of the new model is (RixRj). 

Application of various embodiments of methods and systems to representative, large- 
15 scaled CFD models will now be described. Unlike the general system described by 
Equations (1) and (2), an unsteady fluid dynamic system is driven by displacement and 
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velocity of its moving boundary surface simultaneously as they both contribute to the normal 
downwash on the surface. If one considers a statically nonlinear, dynamically linearized flow 
field, the unsteady fluid motion can be described as 



a? 14 * = Ax* 1 + B Q u^ + Eaii* - (58) . 

y 8 = «(Cx? + D n i^ + Diti l ) (59) 

x = (X x 2) fixM states (60) 
jig s (A^x 1) g^uerafiaed disp&ic^meaifcs (61) 
m s (A^ X 1) geaijera&ad! vdbciies (62) 

y = (Nz X 1) genjeral^di aerodpiajiijc forces (63) 
q = . dominie pressure (64) 

It is noted that the above equations progress in non-dimensional time, r ~~ fe , rather 
than in the real time t and ( ) is the first derivative with respect to t. In fact, the dependency 

10 of the normal downwash on air speed disappears when the equations are discretized in t, as 
in Equations (58) and (59). The structural degrees of freedom, w, are the generalized 
coordinates associated with structural modes. These modes are used to describe the motion 
of the lifting surface. Two different types of reduced-order fluid dynamic models can be 
obtained depending on how the inputs are treated. If necessary, the FDKL/SCI can be 

1 5 performed for additional reduction. 

In one particular embodiment, a method in accordance with the invention may be 
applied to an aerodynamic ROM with displacement and velocity inputs. In this embodiment, 
one can treat w n and u n separately as independent inputs. Thus, for the pulse inputs 
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for i = 1, 2,. . .,7V,, we obtain the corresponding responses rfsJi at the first two time 
steps. Let us define 

— |y{ • ■ • yij yifH i taws - ■ • ystal C^ 7 ) 

where the first JVj* samples are due to the pulses in w n and the next Ni ones are due to 
the pulses in u n . Next, we prepare the following inputs, 

h§cr = E^+E^^-w (68) 



Subject to the SCI, we sample the system response y" and get 
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^Vj 2iVj 

jg,. b y"+» - £ y?rf * - £ rfrf (71) 

sWL s-i 



Define 



Hca b ^ 0 £ ... ^"H (72) 

He! = - ^"H (T3) 

where SVD of yaeMs 

- U*sjf 2$ 3 V£ (74) 

wKbi i? s ranfcCH^D)- HieiKDe. tine i^Tradk?rd«x model is given ly 

= Oe first lb oakums of Y* (75) 

Di = the seoauli JVs coSiuihis of Y 11 (76) 

C - U«S^ (77) 

ft, ^ tlififi^jY^aaSvumisafS^USY 1 (78) 

B L the s»seond JV- QcAamis oES^^Uj^Y 1 (79) 

A a S^US^V^S^ (80) 



In yet another particular embodiment, a method in accordance with the invention may 
be applied to an aerodynamic ROM with only the displacements as the system inputs. This is 
5 possible by applying simultaneously the pulse and double pulse inputs, 



f £ C»-0> J 

u?=^ = 4 -i (n = l) J (82) 
I 0 (n- J 
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and get the corresponding responses Jd *~ m Jdi at the first two time steps: 



For a mw SQ. we use 



- lySa *S 2 - jSas! (83) 

= b^a yia --- (84) 



JYi AS 

= EVurf + Efcwfl (85) 

&m\ ami. 



where are discrete-time derivative of . To be consistent with the use of the 
double pulse defined in (82), *F are obtained by filtering T via i.e., 

= cww(rf B dS) (87) 
which is equivalent to the backward difference equation, 

f? = r * Ar r ° (SS) 

Subject to the new SCI we sample the system response y" and get 
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0 0 







(89) 
(90) 

(91) 
(92) 



(93) 



with R=rank(Hdco)- The new reduced-order model has only TV,- input channels and is in the 
form 



x? 11 = A*" + Bu^ (94) 

= q(Ci? + I>TP) (95) 

where 

D = TfS (96) 

C ^ U«S^ 2 (97) 

B ^ E^I^Yi (98) 

A ^ E^I^H^VflS^ (99) 

Embodiments of aeroelastic systems in accordance with the present invention can be 
constructed using the embodiments of reduced-order aerodynamic models described above. 
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For example, in one embodiment, we first note that structural model is normally described in real, 
continuous-time: 



mia + cti + = y (100) . 

5 

( ) and ( " ), respectively, represent the first and the second derivatives with respect to t. 

Hence, to construct aeroelastic model the continuous-time equation (100) is discretized in 
time: 

2^ - A^ + B^y* (HOI) 

u m = G B iF (102) 

z s { I | (H03) 

a = p oj (ao4) 

10 Note that given At and Fthe consistent incremental time step — must be used in the 
conversion to the discrete-time. 



Aeroelastic Model I 

In this approach, we treat the displacement and velocity, w n , u n as independent inputs 
15 and apply the *%cr. d 3cr given by Eqs. (68) and (69) at a reference dynamic pressure, 
*^ " S*WW corresponding samples are taken and scaled by The ROM of the 
first kind described above is then obtained by applying the SCI/ERA. An aeroelastic system 
that is valid at all V can be constructed by combining the resulting aerodynamic ROM with 
the structural equations: 
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X ra41 = AiXf (105) 
? ~ 1 (106) 



-{:} 



Denoting the eigenvalues of system matrix (107) as ^* h f the aeroelastic eigenvalues 
in the Laplace domain are obtained as 

= ^ (108) 

5 For flutter instability, one must have Real(. > 0 and I— 1| ^'\\ > 1 for any i. 

Aeroelastic Model II 

One can also apply the given by Equations (85) and (86) and get the ROM 

of the second kind described above. This will produce aerodynamic system matrices, A, B, 
10 C, and D where B, D each has only JV/ columns. The new resulting reduced-order aeroelastic 
model can be obtained as 



X** 41 = A^X* (109) 



a - [ A BC S I 

^ ~ [qB a G A*+qBJ)C B \ 



(110) 



Although q can change in Equation (110), this model must be used only at the 
15 reference air speed V re f Given a local angle of attack w, plunging rate u , and local air speed 
V the total aerodynamic downwash at the moving boundary is u + V u and as such it is 
impossible to separate out and account for the effect of V without having both the 
displacement and velocity channels. However, this drawback is easily remedied by adjusting 
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the incremental time step according to when the air speed changes from one value 

to another and discretizing the structural model based on the new At. That is, if one leaves the 
V dependency in the structure, 



■[ 



A BQs 
qB s (V) C A.00 + <?Bs(\ ODC, 



(211) 



then the Aeroelastic Model II becomes valid for all air speeds. 

To demonstrate an embodiment of the present invention, an unsteady vortex lattice 
aerodynamic model subjected to several structural mode inputs may be considered. For 

10 example, FIGURE 3 is a schematic view of a vortex lattice formulation 300 for computing an 
unsteady flow field around a flat, rectangular wing 302 in accordance with an embodiment of 
the present invention. In this embodiment, the unsteady flow field is modeled as an 
incompressible, subsonic air flow. In this example, the wing 302 is 3 inches wide and 12 
inches long, has 10 and 20 vortex elements in the chordwise and spanwise directions, 

15 respectively. A free wake 304 has 40 and 20 vortex elements in the streamwise and spanwise 
directions, creating a total of 800 degrees of freedom, (see Kim, T., Nam, C, and Kim, Y., 
1997, Reduced-Order Aeroservoelastic Model with an Unsteady Aerodynamic Eigen 
Formulation, AIAA Journal, Vol. 35, No. 6, pp. 1087-1088, incorporated herein by 
reference). The wing structure is modeled using 6 vibrational (3 bending and 3 torsional) 

20 modes (see Crawley, E.F., and Dugundji, J., 1980, Frequency Determination and Non- 
Dimensionalization for Composite Cantilever Plates, Journal Sound and Vibration, Vol. 72, 
No. 1, pp. 1-10, incorporated herein by reference). 

No structural damping was introduced in this example. Thus, the size of the full-order 
aeroelastic model is (812 x 812). 
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For the Aeroelastic Model I the reference air density and speed 306 were set at 

I. 23kg/m\ 80m/sec, respectively. The incremental time at this reference speed is ^ "" x ™* 
9.525 x 10 s sec. For the sampling of the vortex model, 480 extra outputs were extracted in 
addition to the 6 generalized aerodynamic forces at 481 time steps. Applying 12 sets of 

5 random signal inputs simultaneously, 6 for w n , 6 for u n , yielded a single set of sampled data. 
12 pulse inputs were also applied individually at the first two time steps to generate Y° and 
Y 1 . 

FIGURE 4 shows three sets of random numbers 402, 404, 406 generated as inputs to 
the vortex lattice formulation of FIGURE 3 using the well-known MATLAB computer 

10 program. Out of 481 time samples, the SVD produced 413 linearly independent singular 
modes. This number was determined by the rank of H c0 matrix. Thus, the size of the 
aerodynamic ROM became (413 x 413). The reduced-order aerodynamic model was then 
coupled with the structural model to create (425 x 425) aeroelastic model (ROM I.). 

For the Aeroelastic Model n. the reference air density and speed were again set at 

15 1.23kg/m 3 , and 80m/sec. 6 sets of random signals and 6 sets of discrete-time derivatives of 
the random signals were applied for u n and u n using 481 time steps and the 486 output 
measurements. This yielded (329 x 329) aerodynamic ROM which when combined with the 
structural system, produced (341 x 341) aeroelastic model (ROM II.). It is noted that ROM 

II. is approximately 20% smaller than ROM I. as a result of using only the half of the input 
20 channels. 

Next, the dimensions of the reduced-order aerodynamic models were further de- 
creased using the FDKL/SCI method. As mentioned earlier, the incremental time step 
embedded in both the FOM and SCI/ERA ROM is too small to be effective for various 
aeroelastic simulations which usually involve a low frequency range. Considering that the 
25 highest free vibrational frequency of the structural modes is 4,160 rad/sec, the sampling 
range in the FDKL method was restricted to (-4,500, 4,500) rad/sec. For the ROM I., out of 
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174 frequency samples within the range 129 KL modes were selected based on the rank of 
the covariance matrix K. Hence, the size of the new reduced-order aerodynamic and 
aeroelastic models (ROM I.-FDKL) became 129 and 141, respectively. Likewise, for the 
ROM II. 97 KL modes out of 130 frequency samples in the same frequency range were 
5 selected yielding new (109 x 109) aeroelastic model (ROM IL-FDKL). For computational 
efficiency, these reduced-order models are to be preferred over the ROM I. and ROM II. 

FIGURE 5 shows three sets of generalized aerodynamic forces 500 computed using 
the vortex lattice formulation of FIGURE 3. Specifically, FIGURE 5 shows (6 x 6) 
generalized aerodynamic forces for V=80 m/sec obtained from the FOM, ROM I.-FDKL, and 
10 ROM IL-FDKL models (FOM (800), ROM I.-FDKL(129), and ROM IL-FDKL(97)), in the 
nondimensional time, t. It is seen that despite the cut-off frequency range present in the 
latter two models, they reproduce the pulse aerodynamic responses of the original model very 
well. 

FIGURES 6 and 7 show two different scales of the aeroelastic eigenvalues 600, 700 
1 5 of the various models in the s domain at V = 80m/sec. It is seen that many eigenvalues of the 
reduced-order aeroelastic models match very well with those of the full model (FIGURE 6). 
More specifically, the 12 complex eigenvalues associated with 6 structural modes agree very 
well between the FOM and ROMs, although the higher structural modes (5th and 6th) in the 
ROM H. and ROM IL-FDKL are slightly mismatched (FIGURE 7). Also noteworthy is that 
20 all the eigenvalues of the ROM I.-FDKL and ROM IL-FDKL are approximately within the 
specified bound, (-4,500, 4,500) rad/sec as the model was obtained based on frequency 
samples in the range. 

FIGURES 8 and 9 show time responses of the first two structural modes 800, 900 due 
to an initial condition in tti. It can be seen that the three sets of curves are practically 
25 indistinguishable from each other. 
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Similarly, FIGURES 10-13 show the aeroelastic results at V = 121.2 m/sec. 
Specifically, FIGURE 10 shows a set of aeroelastic eigenvalues 1000 computed using the 
vortex lattice formulation of FIGURE 3 at a first scale for V = 121.2 m/sec. FIGURE 11 
shows a set of aeroelastic eigenvalues 1100 computed using the vortex lattice formulation of 
5 FIGURE 3 at a second scale for V = 121.2 m/sec. FIGURE 12 shows a time response 1200 
of the first structural mode computed using the vortex lattice formulation of FIGURE 3 for V 
= 121.2 m/sec. And FIGURE 13 shows a time response 1300 of the second structural mode 
computed using the vortex lattice formulation of FIGURE 3 for V = 121.2 m/sec. As can be 
seen from these figures, the wing is on the verge of flutter at this speed. It may be noted how 

10 accurately the ROM I.-FDKL is able to reproduce neutrally stable, sinusoidal time responses 
(FIGURES 12-13). However, the ROM II.-FDKL exhibits a noticeable but minor error in 
producing the transient response. 

FIGURE 14 shows a graph 1400 of model construction time of VLM ERA ROMs 
versus number of inputs for the vortex lattice formulation of FIGURE 3. Using FIGURE 14, 

15 the model construction time may be compared between the two ERA methods. In order to 
obtain accurate and consistent singular modes, H c0 , H c i matrices were kept as square as 
possible by keeping the number of time samples approximately equal to the total number of 
measurements which is the sum of the number of generalized aerodynamic forces and the 
number of auxiliary measurements. The same numbers of time steps and auxiliary outputs 

20 were used both in the Pulse/ERA and Sd/ERA. Thus, in the first case where only the first 
bending mode alone excites the flow field, M = 131, N 0 = 131. In the second case where the 
first bending and first torsional modes were included, M = 251, No = 252, and in the third 
case where the first bending and torsional as well as the second bending modes excite the 
aerodynamic field, M = 281, N 0 = 283. 4 and 5 inputs were also used with M = 331, N 0 = 334 

25 and M = 411,Afo = 415, respectively. 
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FIGURE 15 shows CPU seconds spent constructing various models in accordance 
with embodiments of the present invention (Note that FIGURE 14 was obtained based on the 
data in FIGURE 15). Specifically, Table 1 of FIGURE 15 shows CPU seconds spent in 
constructing ROM I. on a SGI machine. Also presented in parenthesis are the dimensions of 
5 the corresponding reduced-order models. Table 2. shows CPU seconds consumed for ROM 
II. on the SGI machine. The numbers shown in FIGURES 14 and 15 represent total CPU 
seconds spent not only in sampling the response but also processing the data in the 
subsequent ERA schemes. As seen in FIGURES 14 and 15, the new method clearly has an 
advantage over the Pulse/ERA in reducing the model construction time yielding saving 

10 factors of multiple numbers. Needless to say, as the number of inputs increases so does the 
saving. It is interesting that for a given number of inputs both ERA methods generate ROMs 
of very similar sizes. As expected, ROM II of the SCI/ERA are as small as 80 % of the 
corresponding ROM I. Despite the different input channels, both SCI/ERA ROM I. and 
ROM n. require approximately the same CPU time implying that in the SCI/ERA the overall 

1 5 computation time is not very sensitive to the number of inputs. 

In yet another embodiment of the present invention, the SCI/ERA method has also 
been applied to unsteady aerodynamic systems modeled by CFL3D code. CFL3D is a finite 
element program that based on the Navier-Stokes equations models nonlinear viscous, 
compressible fluid motion in subsonic as well as transonic flow fields, (see Krist, S.L., 

20 Biedron, R.T., and Rumsey, C.L., CFLSD User's Manual (Version 5.0), The NASA Langley 
Research Center, Hampton, VA). Although CFL3D describes statically and dynamically 
nonlinear flow, when subjected to a small amplitude it predicts dynamically linearized 
behavior of the flow field around the nonlinear static position in the form of Equations (58) 
and (59). For most practical aeroelastic analyses such as flutter prediction and dynamic gust 

25 loads calculation, the small amplitude approximation yields sufficiently accurate results. 
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The aeroelastic formulation for the CFD code is slightly different than for the vortex 
lattice case in that the second aerodynamic model discussed above is used but whenever the 
air speed V is changed the incremental time step At is adjusted accordingly when discretizing 
the structural model, Equation (100). As in the Aeroelastic Model I, the resulting model can 
5 account for the effect of changing the free stream speed. 

For example, FIGURE 16 is a representative CFD simulation model 1600 in 
accordance with an embodiment of the present invention. In this embodiment, the airplane 
configuration under consideration is that of the Twin-Engine-Transport-Flutter-Model 
(TETFM). The aerodynamic grid is given by the so called the Wing-Pencil-Nacelle (WPN) 
10 model with the strut between the wing and nacelle omitted. The structural motion is 
described by 10 structural modes, resulting in total of 10 generalized aerodynamic forces per 
a mode shape. 

The computational aerodynamic model consists of approximately 700,000 cells and 
30 blocks. For detailed description of the modeling, see Hong, M.S., Bhatia, K.G., SenGupta, 

15 G., Kim, T., Kuruvila, G., Silva, W.A., Bartels, and R., Biedron, R., Simulations of a Twin- 
Engine Transport Flutter Model In the Transonic Dynamics Tunnel, IFASD Paper 2003-US- 
44, incorporated herein by reference. Also, the details of the computational model and 
construction of the aerodynamic and aeroelastic ROMs based on the SCI/ERA using various 
types of input signals mentioned earlier is described more fully in Kim, T., Hong, M.S., 

20 Bhatia, K.G., SenGupta, G., Aeroelastic Model Reduction for an Affordable CFD Based 
Flutter Analysis, AIAA Paper 2004-2040, published subsequent to the filing of the present 
application and incorporated herein by reference. 

FIGURE 17 shows four structural mode shapes 1700 used for the CFD simulation 
model 1600 of FIGURE 16. 

25 FIGURE 18 shows a V-g plot 1800 of the two different aeroelastic models in 

accordance with the previous and an embodiment of the present invention, respectively. The 
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WPN model was examined at M = .831. With A/ = 3.34 x I0~ 4 sec and 995 time steps, 
aerodynamic ROMs were created using both the Pulse/ERA and SSCI/ERA methods based 
on the displacement method described earlier (ROM II.). The size of the aerodynamic ROMs 
is (512 x 512). Given the 10 mode inputs, the total number of time steps consumed in the 
5 time marching required for the SGI/ERA process was 995 + 81 x 10 = 1805. On the other 
hand, the traditional ERA required total of 10 x 995 = 9, 950. Using a single CPU of 
IBM/Regatta machine, the total CPU hours for the Pulse/ERA calculation was 366 while the 
SSCI/ERA required only 64 hours. Thus, the computational cost was reduced by a factor of 
5.7. As shown in FIGURE 18, the two ROMs match very well in flutter characteristics 

10 predicting the two flutter points within just 0.13 % and 0.36 % differences, respectively. 

In accordance with embodiments of the present invention, efficient time-domain 
model reduction/system identification techniques have been presented and demonstrated for 
linear dynamic systems that are subjected to multiple right-hand-side inputs. Methods and 
systems in accordance with the present invention do not require modifying the original code 

15 and take only input and output data for the model construction. Furthermore, methods in 
accordance with the present invention are based on a direct singular value decomposition of 
the output measurements that are not necessarily attributed to pulse inputs but due to multiple 
signal inputs applied simultaneously at the input channels. Compared to the Pulse/ERA, the 
SCI/ERA methods disclosed herein can significantly reduce the model construction time and 

20 compress the amount of output data. Therefore, such methods and systems are very attractive 
for large-scaled dynamic systems with multiple driving inputs such as CFD models wherein 
the moving boundary input is often described by many structural modes. 

While preferred and alternate embodiments of the invention have been illustrated and 
described, as noted above, many changes can be made without departing from the spirit and 

25 scope of the invention. Accordingly, the scope of the invention is not limited by the 
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disclosure of the preferred and alternate embodiments. Instead, the invention should be 
determined entirely by reference to the claims that follow. 



.36. Black Lowe & Graham F 

25315 ^ 

CUSTOMER NUMBER BING-l-lWOAP.doc 701 Fifth Avenue, Suite 4800 

Seattle, Washington 98104 
206.381.3300 • F: 206.381.3301 



